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In the overall aim for a better understanding of the 
vestibular and optokinetic systems and their roles in space 
motion sickness, the eye movement responses to various dynamic 
stimuli are measured. The vest ibulo— ocular reflex (VOR) and 
the optokinetic response, as the eye movement responses are 
known, consist of slow phase and fast phase nystagmus. 

The specific objective of this study is to develop 
software programs necessary to characterize the vestibulo- 
ocular and optokinetic responses by distinguishing between the 
two phases of nystagmus. The overall program is to handle 
large volumes of highly variable data (nystagmus waveforms) 
with minimum operator interaction. The programs include 
digital filters, differentiation, identification of fast 
phases, and reconstruction of the slow phase with a least 
squares fit such that sinusoidal or psuedorandom data may be 
processed with accurate results. The resultant waveform, slow 
phase velocity eye movements, serves as input data to the 
spectral analysis programs previously developed for NASA, 

Johnson Space Center, Neurophysiology Laboratory to analyze 
nystagmus responses to psuedorandom angular velocity inputs. 
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Introduct ion 


Stimulation of the vestibular system by angular 
acceleration during head movements results in a reflexive 
eye movement called nystagmus. The resulting response 
(nystagmus) resembles a sawtooth waveform of which slow 
rotation of the eyes to maintain gaze on an object are 
related to the simulus (slow phase) while rapid resets of 
the eye position are related to some centering mechanism 
(fast phase). Quantitative properties of nystagmus are 
important in the characterization of the vestibular and 
ocular systems. Thus* the specific objective of this study 
is to develop the necessary software programs to 
characterize the vestibular and optokinetic responses to 
sinusoidal stimuli by distinguishing between the two 
components of nystagmus. 

Background 

The methods used to analyze nystagmus vary from manual * 
to real-time automated processing in the time domain. 
Massoumnia presented models of the slow and fast phase 
velocities to establish that the slow and fast phase 
velocity spectra were superimposed. Thereby* he concluded 
that it was not possible to distinguish between the slow 
phase velocity and the fast phase velocity by analysis in 
the frequency domain (7). 



In general, there are two time domain methods for 
identification of a fast phase event. One method is based 
on the analytical geometric property that points at which 
the first derivative waveform equals zero correspond to 
point at which the position waveform is an extremum (maximum 
or minimum). This method involves digital filtering and 
differentiating the ocular position signal prior to 
detecting the fast event with some velocity criterium. The 
second method is based on identification of the position 
waveform extrema by a search algorithm. Then a psuedo 
position waveform is obtained by connecting straight line 
segments between maxima and minima. The psuedo position 
waveform is used to calculate slopes (velocities), position 
changes, and time durations which in turn, are used to 
identify the fast phase events (8) . An assumption of the 
method is that the information of value which is contained 
in the nystagmus can be obtained by replacing the detailed 
time-varying path of the EOS with straight-line 
approximations between points of maxima and minima. Seven 
parameter can be extracted by the method as is shown in 
figure 1. To date, one major vestibular research laboratory 
in the United States, out of the seven major laboratories 
which were selected for this study, uses the psuedo waveform 
method . 
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Upper trace shows the psuedo position EOG waveform 
superimposed on the raw data. Lower drawing is 
an expanded version of the straight line approx- 
imation between points of extrema (MAX and MIN). 

SA denotes slow phase amplitude, FA is fast phase 
amplitude, SV is slow phase velocity, FV is fast 
phase velocity, SD is the slow phase duration, 

FD is the fast phase duration and T is the period; 
where frequency is 1/T. 
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Char act eristic Parameters 


In characterizing the slow phase velocity component 
vestibular researchers have used the engineer’s approach to 
the analysis of a system, i.e., observe the response to a 
predetermined excitation in order to determine the system s 
transfer function. Knowledge of the system’s transfer 
function permits prediction of the system response to other 
deterministic excitations. The most commonly used 
excitation signals in signal analysis are single frequency 
sinusoid, step, and impulse functions. The most popular 
analytical method is the use of a single frequency 
sinusoidal stimulus which applies only to linear systems 
analysis. The only parameters that can vary are the 
magnitude and phase of the resulting sinusoid. The 
frequency of the input and output are unchanged unless the 
system is nonlinear. Thus, gain (the ratio of output 
magnitude to input magnitude) and phase (the delay of signal 
transmission through the system) at a constant frequency are 
the only two parameters necessary to characterize the slow 
phase velocity systems. 

Characterization of the fast phase velocity components 
of nystagmus (saccades) is not obtained by the frequency 
response method but rather by identifying true saccades, 
then measuring the maximum velocity of the fast phase 
segment (saccade), the total change of amplitude from the 
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start to the end of the fast phase segment, and the duration 
of the saccade. A least square exponential curve fit to the 
saccade maximum velocity versus total change of the saccade 
amplitude data results in two best fit coefficients. In 
addition saccade velocity, accuracy, and reaction time were 
used test the oculomotor system (2, 3, A and 5). 

Current Processing 

A general review of the analog and digitizing processes 
used by the laboratory reveals no agreement on filtering or 
digitizing as seen in Table I. Most laboratories used DEC 
LSI-11 computers with 12 bit analog-to-digi tal converters. 
The sampling rate at which the EQG data is digitized depends 
upon the objective of the analysis. If the principle goal 
is to identify and quantify the characteristics of the fast 
phase velocity components of nystagmus (saccades) as a means 
of diagnosing vest ibulo- ocular disorders (3) , then the data 
is sampled at 200 samples per second. Whereas, if the 
primary interest is solely on removal of the fast phases as 
a means of obtaining a more accurate estimate of the slow 
phase velocity parameters, then lower sampling rates are 
used . 

The algorithms for processing the eye movements vary in 
structure primarily because of the ultimate aim of the 
analysis. Laboratories interested in retaining fast phase 
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TABLE I 


Comparison of Analog Antialiasing 
Filtering and Digitizing 
in Vestibular Laboratories 

LAB 


ANALOG 

A 

B 

C 

D 

E 

F 

G 

* 

PREAMPS 

IH 

IH 

IH 

IH 

IH 

IH 

IH 

FILTER 

Bu 

6-pt 

2-pt, 

none 

Bu 

Bu 


TYPE 


Bessel 

Bu 





CUTOFF (Hz) 

35 

41.6 

80 

— 

35 

35 


GAIN 

— 

— 

— 

— 

— 

— 

— 

DIGITAL 








COMPUTER 

PDP 

LSI 

LSI 

LSI 

LSI 

LSI 

LSI 


11/34 

11 

11 

11/73 

11/73 

11/73 

11 

A/D (Bits) 

12 

12 

12 

12 

12 

12 

12 

SAMPLING RATE (S/S) 

122. £ 

200 

200 

200 

200 

100 

120 


IH Fabricated in-house 
Bu Butterworth filter 
S/S Samples per second 
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information use higher bandwidth digital filters in order to 
maintain waveform and timing accuracy. 

Two laboratories use optimal band limited derivatives 
(BLD) in lieu of some smoothing routine followed by a two- 
point central difference equation (CDE). The two-point 
central difference is popular as a first order 
differentiator because of its speed, simplicity, accuracy 
and low pass filtering (1). The laboratories using optimal 
filtering techniques convolve a finite impulse response 
(FIR) filter with an finite impulse differentiator to obtain 
velocity or acceleration filters. TABLE II compares the 
digital processes for eye movement analysis used by the 
seven vestibular laboratories. 

It is interesting to note that three laboratories 
identify three classes of fast phase events but only one 
laboratory uses the saccade velocity information. In the 
evaluation of the slow phase velocity, most laboratories 
have algorithms which remove the fast phase events and fill- 
in the removed points with a linear extrapolation across 
either the position or velocity waveform. Two laboratories 
use a least squares sinusoidal fit to the velocity curve 
without filling the gap between slow velocity segments. 

Three laboratories have equally spaces intervals after the 
waveform reconstruction and only two of those laboratories 
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TABLE II 


Comparison of Digital Processes 
For Eye Movement Analysis 
Laboratory 



A 

B 

c 

D 

E 

F 

G 

FILTER TYPE 
CUTOFF (Hz) 

BLD 

None 

7-pt 

LP 

20 

LP 

ReF 

35 

BLD 

None 

15-pt 
LP FIR 
25 

DERIVATIVES 





* 



1st 

2nd 

— 

2-pt 

CDE 

2-pt 

CDE 

3-pt 

CDE 

3-pt 

CDE 

2-pt 

FWD 

4-pt 

CDE 

VEL. FILTER 
CUTOFF 

31-pt 

25 

— 

1 

— 

9-pt 

FIR 

— 

— 

ACCEL. FILTER 
CUTOFF 

61-pt 
3Q Hz 

— 

— 

— 

— 

— 

— 

EXTRAPOLATION 

Lin 

Lin 

None 

None 

Lin 

Lin 

Lin 

SLOW PHASE WAVE 
RECONSTRUCTION 

Vel 

Vel 

None 

None 

Vel 

Pos 

Pos 

FAST PHASE 
CLASSIFICATION 


3 

— 

3 

3 

— 

— 

LEAST SQUARES 
TRANSFORM 

— 

— 

Sin 

Avg 

SPV 

— 

Sin 

-- 

DIRECT FOURIER 
TRANSFORM 

— 

yes 

— 

— 

— 

— 

— 

FFT 

j yes 

— 

— 

— 

— 

— 

yes 


CDE 

Central difference 

equation 

Sin 

Sinusoidal 

FWD 

Forward difference 

equation 

Pos 

Position 

FFT 

Fast Fourier Transform 

Lin 

Linear 

Vel 

Velocity 


ReF 

Recursive Filter 
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use a fast fourier transform to obtain the frequency 
response parameters. 

NASA FPID Program 

Development of the fast phase identification <FPID) 
program began prior to evaluation of the seven U.S. 
laborator ies . The engineering approach dictated the 
analytical geometric method. Horizontal eye movements (VOR) 
response to a sinusoidal stimulation of the vestibular 
system is digitized at 120 samples per second (Fig 2) and 
filtered with a digital 15-point, low-pass, finite-impulse- 
response (FIR) filter. The FIR filter cutoff is set at 25 
Hz. The signal at 36 Hz is -^0.1 db . Several FIR filters 
and smoothing routines were evaluated before selection of 
the final 15-point, FIR filter. The filtered EGG (position) 
signal is shown in figure 3. 

The filtered signal is then differentiated with a 
central difference equation. The program permits the 
operator to select the polynomial order of the 
differentiator up to a sixth order central difference 
equation with error of order (h A ), where h is the sampling 
interval of 0.00833 seconds (8.33 msec.). Higher order 
equations result in better accuracy, but increase 
computation time. Hence, the fourth order <A-point) central 
difference equation with error of order (h“*) is used to 
obtain the first derivative of the EOG since it gives 
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shown after 



FIR FILTERED DATA 
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Fig. 3. Horizontal eye movements after FIR filtering. 


accurate results with only two points on either side of 
center (fig. *f) . 

The second derivative is obtained by differentiating 
the filtered EOG position signal with a 5-point, second 
derivative central difference equation of error order (h^). 
The second derivative is shown in figure 5. The third 
derivative of the EOG position signal was also obtained with 
a 7-point central difference equation (6). The third 
derivative waveform was noisier than the lower derivative 
filters so it was removed from the program. 

Following differentiation, the program computes the 
root-mean-square (RMS) value of the various derivatives 
waveforms in order to set a threshold value above the noise 
level so as to reduce false detection errors. From RMS 
values and approximate signal-to-noise ratios of each 
derivative waveform, detection of a fast velocity event is 
based on exceeding the first derivative threshold value 
rather than the second derivative as used by Massoumnia (7). 
The first derivative offersd the least amount of noise and 
the best signal— to— noise ratio for threshold detection. 

Once the threshold of the EOG velocity signal is 
exceeded the search point is move backward to find where the 
derivative zero crossing occurrs. This point is flagged as 
the beginning of the fast phase event and the search is 
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FIRST DERIVATIVE OF DATA 



First derivative of the filtered eye mov 
The RMS value is 585 A/D counts for 2.08 
of data. The signal-to-noise ratio is 3 




SECOND DERIVATIVE OF DATA 
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reversed (forward direction) until a zero crossing occurrs. 
This point was flagged as the end of the fast phase event. 

In the next steps, the program uses the filtered eye 
position (EQG) waveform to perform a least squares linear 
regression on the slow phase velocity segment preceding the 
starting index of the current fast phase event. Then the 
points between start and end of the fast phase event are 
extrapotated and added to the position waveform. 

Prior to output of the reconstructed slow phase EQG 
signal, a correct ion * of slow segment height is necessary. 
This is accomplished by obtaining the change of position 
(height) from the last point of the extrapolated EGG 
position (end of fast phase index) and the start of the 
following slow phase segment (end point + 1). Height 
corrections are cummulative and must carry the proper sign. 
The r econs tr uc t ed slow phase EOG position waveform is shown 
in figure 6. The FPID program listing is given in the 
appendix . 

Conclusion and Recommendations 

Although the fast phase identification (FPID) and slow 
phase reconstruction appear to work, the program needs to be 
evaluated with various types of VOR and OKN data. At the 
present time only short segments of data can be analyzed. 

The program needs a circular buffer, so that long data 
records can be analyzed. Additionally, the reconstructed 
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RECONSTRUCTED SLOU PHASE EOO POSITION DATA 



Fig. 6. The reconstructed slow phase horizontal eye 


wave must be written to a new data file prior to being used 
by the spectral programs written for NASA last summer. 
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APPENDIX 


FAST PHASE IDENTIFICATION 
( FPID ) 

PROGRAM LISTING 
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ORIGINAL PAGE IS 
OF POOR QUALITY 


F u K f Kf-iN I l J VO 2 * 6 Weo 06-tfua-S6 1 7 ; 0 X : t>4 


PrtGE O0L 


0001 

C 

0002 

0003 

C 

0004 


0005 

0006 

0007 

0008 

0009 

0010 
0011 
0012 

0013 

0014 

0015 

0016 5 
C 

0017 

0018 

0019 

c 

0020 
0021 

0022 10 

0023 20 

0024 

0025 30 

0026 

0027 40 
C 

0028 
0029 

C 

0031 42 

0032 

0033 45 

0034 

C 

C 

0035 

0036 

0037 

C 

C 

0038 

C 

0039 

0040 


F'RuGftrtfl hP ID 

INTEGER Z( 512) . ASKING > NBLK » IPOSN ( l») ,MAX< lo) . ICHNUM< lo / .LENGTH 
REAL AMPL. PERIOD. INTRaT. HF0G<512). XT',512). COEF(IS). 

1 YHE0G(512). SDE0G(512>» SPEOG ( 512 .< . SP YHFOG ( 512 > . FOEOGi5l2) 

EQUIVALENCE ( AriPL » Z ( 150) ) » 1 PERIOD. Z< 152) ) » i I nTRaT. H 15?) ) , 

1 ( DSKINC » Z < 167 ) ) » (NBLK . Z < 172 ) ) » ( ICHNL'M . Z ( 1 77 ) ) . 

2 (IPOSN. Z(225)). (MAX. Z<241>) 

COEF(l) = .0037165603 

COEF ( 2 ) = .020235427 

COEF < 3 ) = . 013399956 

COEF ( 4 ) = -.040643737 

COEF ( 5 ) = -.066073574 

COEF ( 6 ) = .061152663 

COEF ( 7 ) = .30294424 

COEF 8 ) = .43047235 

NEH = 8 

DO 5 N = 1 , 7 

COEF < 16-N) = COEF < N ) 

CONTINUE 

HHT = 0. 

THT = 0. 

ISLG = 0 

INTEGER 6UFR2 ( 1024 ) . IERR 
BYTE FILNAM ( 12 ) 

DO 20 N*1 » 12 
FILNAM ( K ) =0 
TYPE 30 

FORMAT (' ENTER COMPLETE FILENAME'/' FILENAME ?'.$> 

ACCEPT 40. FILNAM 
FORMAT (12A1) 

CALL D ISNIO ( FILNAM . -3 . Z . 256 . 0 . NDUMMY » IERR ) 

IF (IERR .NE. 0) TYPE *.' ERROR CODE'. IERR.' DURING READING' 

TYPE *» ' CHAN NO. ' 

ACCEPT *.K 

TYPE *»' ENTER THE NUMBER OF DATA POINTS FOR EOG ' 

ACCEPT *. LENGTH 

TYPE *.' DO YOU WANT FILES LEFT OPEN? <1*YES)' 

ACCEPT *.LOPEN 
LOPEN * 1 

TYPE ENTER STARTING BLOCK NUMBER' 

ACCEPT #.IBLN 

TYPE DO YOU WANT TO PLOT EVERY POINT? ENTER STEP SIZE!' 

ACCEPT *.ISTEP 
ISTEP = 1 

NURDS = DSKINC* 256 
I TEMP = 0 
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FORTRAN IV 

Y02.6 Wed 06 -Aijc>-S<j l7J0i:54 

0041 

43 

TYPE 50, NBLK, NURDS* IBLK » DSKINC 

0042 

50 

FORMA T (' NBLK* ' » 15 » ' NURDS* ' » 1 5 » ' IBLK = ' , 15 , ' DSKINC 

0043 


TYPE 51 .MAX IK > f IPOSN(K) 

0044 

51 

(T* 

FORMAT (' MAX(K> * ',15,' IPOSNCK) = '»I5> 

0045 


IRmCDE*3 

0046 


IF (LOPEN .EG* 1) I RHODE * -3 

0048 

c 

TYPE *f' NOW READING FROM DISK ' 

0049 

52 

CALL DISKIO(FILNAM, I RHODE , BUFR2 , NURDS » I BLK , NDUMM it I - j 

0050 

C 

IF < IERR . NE . 0) TYPE ERROR CODE' » I ERR, ' DURING~f 

0052 


NPTS = (MAX<K)-IP0SN(K)+1)/I3TEP 

0053 


DO 57 L = If NPTS 

0054 


HEOG< L+ITEMP ) * BUFR2<IP0SN<K) + L*I5TEP> 

0055 

57 

C 

CONTINUE 

0056 

DO 950 1=1.26 

0057 


TYPE 951. I. HEOG < I ) 

0 058 

951 

FORMAT ( ' I = ' , 14 » ' HEOG = '.F12.5) 

0059 

950 

CONTINUE 

0060 


DO 960 I = LENGTH-26. LENGTH 

0061 


TYPE 961. I. HEOG(I) 

0062 

961 

FORMAT ('I * ' » 14 » ' HEOG = ',F 12.5) 

0063 

960 

CONTINUE 

0064 


TYPE 962. NPTS. I TEMP 

0065 

962 

FORMAT (' NPTS * '.18.' ITEMP = ',13) 

0066 


TYPE WAITING TYPE 1' 

0067 

r 

ACCEPT*. ISN 

0068 

w 

IBLK » IBLK + DSKINC 

0069 


ITEMP * ITEMP + NPTS 

0070 

L 

IF (ITEMP .LT. LENGTH) GO TO 52 

0072 


XLGT * FLOAT < LENGTH ) 

0073 

c 

TYPE 50, NBLK, NURDS, IBLK, DSKINC 

0074 

DO 60 J * 1, LENGTH 

0075 


XT(J) = FLOAT ( J) 

0076 

60 

CONTINUE 

0077 


TYPE *, ' SET YMIN ! ' 

0078 


ACCEPT * » YMIN 

0079 


TYPE SET YMAX ! ' 

0080 

p 

ACCEPT * , YMAX 

0081 

L» 

ILFT = 600 

0082 


IRIT * 3500 

0083 


IBOT * 1000 

0084 


I TOP * 2500 

0085 


XMIN - 0. 

0086 

c 

XMAX * XLGT 

0037 

TYPE *»' DO YOU WANT TO REMOVE D.C.9 <YES*1)' 

0088 


ACCEPT *, I l E3 
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FORTRAN IV 
v>069 

C 

C 

<■>091 
0092 
(>09 3 

0094 58 

0095 

0096 

0097 

0098 59 
C 
C 

0099 

0100 
0101 
0102 

0103 

0104 

0105 


C 


0106 

0107 

0103 

61 

r* 

0109 

0110 
0111 

L> 

c 

0113 

c 

0114 

0115 

62 

0116 

500 

0117 

0113 

0119 

0120 

C 

0121 

0122 

0123 

510 


C 

0124 

162 

0125 

C 

0126 

0127 

0128 

0129 

0130 
0132 



COW .6 Ui c- d 06 — Aii£t — 86 1 7 i 0 l : 4 POoE 003 

IF (I YES • NE < 1) GO TO 162 

COMPUTE THE MEAN OF THE DATA 

XSUM = 0.0 

DO 58 I ■ L LENGTH 

XSUM = XSUM + HEuG(I) 

CONTINUE 

XMEAN = XSUM/X'LGT 
HO 59 I * If LENGTH 
HEOG(I) = HEOG(I) - XMEAN 
CONTINUE 


CALL TSXCHK 

CALL GRINI T ( 401 4 » 4631 » 1 ) 

CALL CHRSIZ ( 2 ) 

CALL ERASE 

CALL GRID( 10 » 10. ILFT » IRIT , IBOT * I TuP, 97) 

CALL ANOT AT ( 10 » 1 0 , ILFT , IR I T » IBOT , I TOP , XMI N • XMA X , YMI N , f MAX ) 

CALL XYPLOT( XT.HEOG.LENGTH. ILFT. IRIT, IBuT, ITuP.XMlrJ-xnAX , 

1 YMIN.YMAX. 1 ,0) 

CALL MPLOT < ILFT+500, IT0P+300.-1 > 

TYPE 61 

FORMAT ( ' +RAW DATA ') 

TYPE DO YOU WANT A HARD COPY? (YES=1)' 

ACCEPT #,IS 

IF (IS .NE. 1) GO TO 62 
CALL COPY (0) 

DO 500 1=3, LENGTH-2 

YHEOG(I) = . 1 1* ( HEOG ( 1-2 ) +HEOG ( 1+2 ) ) + . 22# < HEOG ( 1-1 ) +HEOG ( 1+1 ) ) 
1 + . 33#HE0G ( I ) 

CONTINUE 

YHEOG(l) = YHEOG < 3 ) 

YHE0G<2> = YHEOG ( 3 ) 

YHEOG(LENGTH) = YHEOG ( LENGTH- 2 > 

YHEOG ( LENGTH- 1 ) = YHEOG (LENGTH) 

DO 510 1=1, LENGTH 
HEOG ( I ) = YHEOG ( I ) 

CONTINUE 

NCO = 15 

TYPE FILTERING DATA WITH 15-POINT FIR FILTER ! ' 

DO 80 I * 2, LENGTH+NCO 
SUM =0.0 
DO 70 J » l.NCO 
L = J 

IF (J .GE. I) GO TO 70 
IF (I-J , GT • LENGTH) GO TO 70 
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FORTRAN IV 

Y02. <$ wed 06-Aus-S6 1?:C1.*S4 page 

0.134 


H = COEF < J > * HEOGxI-j) 

0135 


SUrt = SUrt + H 

0136 

70 

CONTINUE 

0137 


YHEOG(I-l) = SUN 

0138 

80 

1*' 

CONTINUE 

0139 


J2 = 14 

0140 


DO 81 J a 15» LENGTH+J2 

0141 


YHEOG (J-J2) = YHEOG (J) 

0142 

81 

CONTINUE 

01 43 


DO 630 I = 1, 14 

0144 


YHEOG (LENGTH- I -1 )»YHEOG( LENGTH-1 4) 

0145 

630 

C 

CONTINUE 

0146 

CALL TSXCHK 

0147 


CALL GRINIT ( 4014 , 4631 » 1 ) 

0148 


CALL CHRSIZ(3> 

0149 


CALL ERASE 

0150 


CALL GRID(10,10, ILFT » IRIT , IBOT , I TOP, 97 ) 

0151 


L ALL ANOTAT ( 10, 10 » ILFT > IRIT , IBOT , I TOP , xhIN * XhAX » YnlN , rrtrtX 

0152 

p 

CALL XYPLOTtXT, YHEOG, LENGTH, ILFT, IRIT, IBOT, ITOP,XrtIN,xnAX 
1 YMIN, YhAX ,1,0) 

0153 


CALL rtPLOT ( ILFT+500, ITOP+300,-1 ). 

0154 


TYPE 261 

0 1 55 

261 

r 

FORMAT ( ' + FIR FILTERED DATA ') 

0156 


TYPE DO YOU WANT A HARD COPY? (YES=1)' 

0157 


ACCEPT *, IS 

0158 

p 

IF (IS .NE. 1) GO TO 83 

0160 

w 

C 

83 

CALL COPY ( 0 ) 

0161 

H = FLOAT ( ISTEP)/120. 

0162 


L = LENGTH 

0163 


TYPE SET DESIRED HALF ORDER OF DIFFERENTIATOR 1 ' 

0164 


ACCEPT *,NDIF 

0165 

p 

GO TO (84,86,88) , NDIF 

0166 

84 

DO 85 I = 2, LENGTH 

0167 


FDEOG(I> = ( YHEOG( 1+1 ) - YHE06 (1-1) )/(H *2.) 

0163 

85 

CONTINUE 

0169 


FDEOG(l) = FDEOG ( 2 > 

0170 


FDEOG(L) * FDEOG (L-l) 

0171 

C 

86 

GO TO 1000 

0172 

DO 87 I » 3 , L-2 

0173 


FDEOG ( I ) » ( ( YHEOG ( I -2 ) - YHEOG( 1+2) > + (8.#( YHFOG( 1+1 ) - 
1 YHEOG (I-1))))/(H#12«) 

0174 

87 

CONTINUE 

0175 


FDEOG ( 1 ) = FDEOG ( 3 ) 

0176 


FDEOG( 2) « FDEOG ( 3 ) 

0177 


FDEOG (L) * FDEOG (L-2) 
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vl78 F£mEQGvL~1 > - FDEOG (. L-2 ) 

C 

0179 1000 DO 200 1*3. L-2 

0180 SDEOG < I ) * (( —1 •)*( YHEOG ( 1—2 ) + yHEOGi It2> r < lo • * '. YHEOG ( I + 1 r 
1 YHEOG(I-l))) - 30 . *rh£0G< I ) ;/ t 12 . *n**2 . > 

0181 200 CONTINUE 

0182 SDEOG(l) = SDEOG ( 3 ) 

0163 SDEOG ( 2 > = SDE0G<3) 

0184 SDEOG ( L ) * SDEOG (L-2) 

0185 SDEOG(L-l) = SDEOG < L-2 ) 

C 

0186 REAL XMFSD 

0187 TYPE *t' SETTING MULTIPLYING FACTOR FOR 2ND DERIVATIVE' 

C ACCEPT *. XriFSD 

0188 XriFSD = .01 
C 

0189 DO 400 I » If L 

0190 SDEOG(I) = XMFSD* SDEOG ( I ) 

0191 400 CONTINUE 

C 

0192 FDRriS = 0. 

0193 SSFD = 0. 

0194 DO 700 I = If LENGTH-20 

0195 SSFD * SSFD + FDE0G(I)**2. 

0196 700 CONTINUE 

0197 SSMFD = SSFD/FLOAT (LENGTrt-20) 

0198 FDRriS = SORT (SSMFD) 

0199 TYPE 725 f FDR riS 

0200 725 FORMAT (' RMS OF 1ST DERV. *'fF8.3> 

C 

0201 SDRMS = 0. 

0202 SSSD = 0. 

0203 DO 750 I = 1 f LENGTH-20 

0204 SSSD = SSSD + SDE0G(I)**2. 

0205 750 CONTINUE 

0206 SSMSD = SSSD/FLOAT < LENGTH-20 ) 

0207 SDRMS = SORT (SSMSD) 

0208 TYPE 775 f SDRMS 

0209 775 FORMAT (' RMS OF 2ND DERV. *',F8.3> 

C 

0210 TYPE #f' UAITING ! TYPE 1 TO CONTINUE.' 

0211 ACCEPT* f JES 

0212 GO TO 90 

C 

0213 88 DO 89 I = 4f L-3 

0214 FDEOG(I) = <<YHE0G(I+3) - YHE0G(I-3)> + (9.*( YHEOGt I- 
1 YHEOG( 1+2) ) ) + ( 45. *( YHEOG ( 1+1 ) - YHEOG (I-1))))/(H*6 

0215 89 CONTINUE 

0216 FDEOG(l) * FDE0G<4> 

0217 FDE0G(2) * FDEOG ( 4 ) 

0218 FDE0G(3) » FDEOG (4) 

0219 FDEOG(L) * FDEOG ( L-3 ) 

0220 FDEOG (L-l) » FDEOG (L-3) 

0221 FDE0G(L-2) = FDE0G(L-3> 
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H OK I 

0222 

0223 
022 4 
0223 
0226 
0227 
022 $ 
022 ? 
02^0 
0231 


0232 

0233 

0234 

0235 

0236 

0237 

0239 

0240 

0241 

0242 

0244 

0245 

0246 

0247 

0243 

0249 

0250 

0251 

0252 

0253 


0254 

0255 

0256 

0257 
0253 
0259 

0261 

0262 

0263 

0264 

0266 

0267 

0268 
0269 


original page re 

OR POOR QUALITY 


'An IV V02 . a u&fi OiS'-rt.is-b* l7;oi:S4 


PrtGE 0 


C 

90 


C 

282 


C 

C 

283 


210 

250 

C 

382 

C 

C 

92 

C 


GO TO 1000 


YMIN * YrilrUlO* 

YMAX = Y H A X*# 10* 

CALL TSXCHK 

CALL GRINI T < 401 4 f 4o3l 1 1 ) 

CALL CHRS IZ ( 3 ) 

CALL ERASE 

CALL GRID ( 1 0 f 10 ? ILFT* IRIT » IB07 f I TOR • 97 ) 

LALL AnO i AT ( 1 0 f 1 0 9 1 LFT t [R I f ? X Du r » I TOP r XM IN f ;< riA x • Y n IN • fnAX ) 
CALL XY PLOT < XT t FDEOG ? LENGTH f I LFT * IP I T » I DuT r IT OP * Ah In • a ha .< 9 
X YhINf YMAXf 1 fO> 

CALL MPLGT ( I LFT +500 f ITOP+300 f -1 ) 

TYPE 282 

FORMAT (' + FIRST DERIVATIVE OF DATA ' ) 

TYPE * 9 ' DO YOU WANT A HARD COPY? <YES=1>' 

ACCEPT *fIS 

IF (IS ♦ NE • 1) GO TO 283 
CALL COPY ( 0 ) 

TYPE DO YOU WANT A SINGLE CURVE PLOT? (YES=1)' 

ACCEPT *» I YES 

IF ( I YES .NE. 1) GO TO 210 

CALL TSXCHK 

CALL GRINIT ( 4014 » 4631 ?1) 

CALL CHRSIZ( 3 ) 

CALL ERASE 

CALL GRID(10rl0fILFT»IRIT >IBOT» I TOP » 97) 

CALL ANOTAT ( 10 > 10 t I LFT rIRIT>IBOT» I TOP t XhIN> XMAX » YMIN* YhAX ) 
GO TO 230 
IBOT = 1250 
ITOP = 2750 

CALL XYPLOT ( XT> SDEOGf LENGTH > I LFT » IRIT » IBOT » ITOP » XRlN p XMAX » 
1 YMIN» YMAX » 1 »0 ) 

.CALL MPLOT ( ILFT+500 r ITOP+300 , -1 ) 

TYPE 382 

FORHAT (' + SECOND DERIVATIVE OF DATA') 

TYPE *>' DO YOU WANT A HARD COPY? (YES=1)' 

ACCEPT *,IS 

IF (IS .NE. 1) GO TO 92 
CALL COPY (0) 

TYPE ANOTHER DIFFERENTIATOR? (YES*1)' 

ACCEPT *»IS 

IF (IS «£Q. 1) GO TO 83 
CALL ERASE 

TYPE *»' IDENTIFICATION OF FAST PHASES' 

FDRHSTH * FDRMS#2. 

TYPE 410» FDRHSTH 
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0270 

410 

0271 

C 

110 

0272 


0273 


0275 

112 


C 

0276 

115 

0277 


0278 

411 


C 

0279 


0281 


0282 


0283 


0285 

117 

0286 


0287 

120 

0283 


0290 

1 22 

0291 


0292 

124 

0293 


0294 

412 


C 

0295 


0296 


0298 

126 

0299 


0300 

128 

0301 


0302 

414 


C 

0303 


0304 


0306 

130 

0307 

134 

0308 


0309 

415 

0310 



C 

0311 

135 

0312 


0313 

416 


C 


C 

0314 

136 

0315 


0316 


0317 


0318 


0319 


0320 


0321 



V02.6 Ujcjo OtS'-rt us-S6 l/': 01:54 I* age CO 7 

FORMAT C ' 1ST DERG. RMS THRESHHOLD ~ ' . F . 3 ) 

LTERM = LENGTH - ISLG 
DO 112 J * If LTERM 

IF (A6S<FDE0GiISLG+J)) .GT. ADS c FDftnSTH ) ) GO TO US 
CONTINUE 

ITLST = J + ISLG 
TYPE 411 , ITLST 

FORMAT (' INDEX TSh EXCD AT I =',I5> 

IF (FD£OG< ITLST) .GE. 0.) 60 TO 120 
TYPE *f' NEGATIVE VALUES START' 

DO 117 L » 1. 10 

IF (FDEOG( ITLST - L) .GE. 0.) 60 TO 124 
CONTINUE 

TYPE *>' POSITIVE VALUES START' 

DO 122 L = 1 > 10 

IF ( FDEOG ( ITLST - L) .LT. 0.) GO TO 128 
CONTINUE 

TYPE *r ' NEGATIVE VALUES END' 

ISTRT = ITLST - L 
TYPE 412 • ISTRT 

FORMAT (' STRT NEG INDEX AT I = ',I4> 

DO 126 K * 1» 30 

IF ( FDEOG( ITLST + K> .GE. 0.) GO TO 135 
CONTINUE 

TYPE *f' POSITIVE VALUES END' 

ISTRT = ITLST - L 
TYPE 414 f ISTRT 

FORMAT (' STRT POS INDEX AT I = 'fI4> 

DO 130 LN * If 30 

IF (FDEOG< ITLST + LN) .LT. 0) GO TO 134 
CONTINUE 

ISTP = ITLST + LN 
TYPE 415 f ISTP 

FORMAT <' END POS INDEX AT I ■ 'fI4> 

GO TO 136 

ISTP = ITLST + K 
TYPE 416 f ISTP 

FORMAT (' END NEG INDEX AT I = 'fI4> 

LEAST SQUARES FIT 
XN = 10 . 

DXSUrt = 0. 

DXSSUM ■ 0. 

DYSUM » 0. 

DXYSUM * 0. 

DO 140 ISL * ISTRT-lOf ISTRT 
FI2 » FLOAT ( ISL ) 

DXSUM = DXSUM + FI2 
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'*'02.^ Wed 0 6 - A ' j a - .'5 6 1/ : 0 1 : s .) 

0322 


DXSSUh = DXSSUM + FI2*fM2 

0323 


DY SUM = DYSUM + YHEOG(ISL) 

0324 


DXYSUM * DXYSUH + FI2*YHE0G< ISL) 

0325 

140 

p 

CONTINUE 

0326 


f'ENOM = tXN*DXSSUM - DXSUn*DXSUM) 

0327 


SLOPE ■ (XN#DXYSUM - DXSUM*D YSUM .» / DENGM 

0328 


TINT * (DYSUM#DXSSUh - OXSUM*OXY SUM ) /DEnOM 

032? 


TYPE 4 1 7 r SLOPEr YINT 

0330 

417 

<*• 

FORMAT (' SLOPE ='»F12.5r' YINT ='»F12.5) 

0331 


DO 142 ILF = ISTRT t ISTP 

0332 


XILF » FLOAT ( ILF ) 

0333 


YHEOG(ILF) = SLOPEKXILF + YINT 

0334 


TYPE *t ILF t YHEOG(ILF) 

0335 

142 

P 

CONTINUE 

0336 

Ur 

TYPE *»' TO CONTINUE TYPE: 1 ' 

0337 


ACCEPT# » IX 

0338 

c 

r 

CALL ERASE 

0339 

TYPE *»' RECONSTRUCT SLOW PHASE' 

£\ T A A 
v 

u 

HHT = HHT -T YHEOG(ISTP) -YHEOG< I3TP+1 ) 

0341 


TYPE 420 » HHT 

0342 

420 

p 

FORMAT < ' HEIGHT CORRECTION *'rF12.5) 

0343 

L» 

TYPE 421 » ISLGr ISTP 

0344 

421 

C 

FORMAT (' STRT SP INDX = r'I4r'STP 3 P INDX 

0345 

DO 145 I » ISLG+1 t ISTP 

0346 


SF'EOG < I ) = YHEOG(I) + THT 

0347 


TYPE *t It SPEOG(I) 

0348 

145 

p 

CONTINUE 

0349 

L 

THT » HHT 

0350 


ISLG » ISTP 

0351 


TYPE 425 1 ISLG 

0352 

425 

FORMAT ( ' NEW SEARCH STARTS AT I = 'rI4) 

0353 


TYPE *»' TO CONTINUE TYPE: 1' 

0354 

r 

ACCEPT *r IXS 

0355 

u 

r 

IF (ISLG .LT. LENGTH) GO TO 110 

0357 

u 

TYPE PLOT RECONSTRUCTED SLOW WAVE? (YES 

0358 


ACCEPT *, I YES 

0359 

p 

IF ( I YES .NE. 1) GO TO 149 

0361 

Lr 

YMIN * YMIN*.l 

0362 

c 

900 

YMAX • YMAX*.l 

0363 

CALL TSXCHN 

0364 


CALL GRINIT(4014> 463 1 r 1 ) 

0365 


CALL CHRSI2 ( 3 ) 
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VO 2.6 
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PAGE 009 


0366 

0367 

0368 

0369 


0370 

0371 

0372 

0373 

0374 

0375 


0377 


0378 

0379 

0380 

0382 

0383 
0334 


CALL ERASE 

CALL OR I D ( 1 0 f 1 0 f I LFT f I A: I T • [ DOT > I TOR' > 97 ) 

CALL ANOTAT (IOfIOfILFT f IRITf 180 T , I TOP f XMIN . XMAX f TMIN f YMAX > 
CALL XYPLOT < XT f SPEC'S f LEnGTHf ILFT f IRITf IBOT f ITOPf XMIN f XrlAX f 
1 YMINfYMAXfI fO) 

c 

CALL MPLGT ( ILFT+500f IT0P+300f-1 ) 

TYPE 882 

882 FORMAT (' + RECONSTRUCTED SLOW PHASE EOG POSITION DATA') 

TYPE *f' DO YOU WANT A HARD COPY? (YES*!)' 

ACCEPT *fIS 

IF (IS .NE. 1) GO TO 149 
C 

CALL COPY < 0) 

C 

C 

C262 TYPE #f ' DO YOU WANT 5-POINT SMOOTHING? ( YES=1 > ' 

C ACCEPT # r I YES 

C IF (IYES .NE. 1) 60 TO 149 

C 

C DO 600 I = 3 f LENGTH-2 

C SPYHEOG ( I ) = . 1 1* ( SPEOG ( 1-2 ) +SPEOG < 1+2 ) ) + . 22# ( 5PE0G < 1-1 ) 

C 1 +SPEOG ( 1 + 1 ) ) + . 33*SPE0G ( I ) 

C600 CONTINUE 

C SPYHEOG(l) * SPYHEOG ( 3 ) 

C SPYHEOG (2) * SPYHEOG ( 3 ) 

C SPYHEOG( LENGTH) = SPYHFOG( LENGTH-2) 

C SPYHE0G(LEN6TH-1 ) = SPYHEOG( LENGTH) 

C 

C375 CALL TSXCHK 

C CALL GRINIT<4014f 4631 f 1 ) 

C CALL CHRSIZ ( 2 ) 

C CALL ERASE 

C CALL GRID<10f10fILFTfIRITfIB0TfIT0Pf97> 

C CALL ANOTAT ( 10 f 10 f ILFT f IRITf IBOTf ITOPf XMINf XMAX fYMINf YMAX) 

C CALL XYPLOT ( XT f SPYHEOG f LENGTH f ILFT f IRIT f IBOT f I TOP f XHIN f XMAX f 

C 1 YMINfYMAXfIfO) 

C 

C CALL MPLOT ( ILFT+500f IT0P+300f-1 ) 

C TYPE 82 

C82 FORMAT ('+ FIR FILTERED AND 5-PT. SMOOTHED DATA') 

C TYPE # f ' DO YOU WANT A HARD COPY? < YES=1 > ' 

C ACCEPT *fIS 

C IF (IS .NE. 1) GO TO 149 

C 

C CALL COPY (0) 

C 

149 TYPE #f ' DO YOU WANT A 15-PT FIR FILTER? (YES=1)' 

ACCEPT *fIYES 

IF (IYES .NE. 1) GO TO 150 
C 

DO 100 I * 2f LENGTH+NCO 

DSUM *0.0 

DO 95 J * If NCO 
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002 . 6 

0365 


IF (J .GE. 1) 

0367 


IF (I-J .GT. L 

0389 


DH = COEF(J) # 

0390 


DSUM = DSUM + 

0391 

95 

CONTINUE 

0392 


SPYHEOG ( I -1 ) * 

0393 

100 

CONTINUE 

0394 


DO 101 M = 15f 

0395 


SPYHEOG (M-J2) 

0396 

101 

i* % 

CONTINUE 

0397 

w 

DO 102 I a If 

0398 


SPYHEOG (LENGTH 

0399 

102 

C 

CONTINUE 

0400 


CALL TSXCHK 

0401 


CALL GRINIT ( 40 

0402 


CALL CHRSIZ ( 3 ) 

0403 


CALL ERASE 

0404 


CALL GRID(10fl' 

0405 


CALL ANOTAT ( 10 

0406 


CALL XYPLOT ( XT 


Neti 0 6 - A u 3 - 6 4 1 7 : 0 1 S4 

130 TO 93 
EN6TH > GO TO 93 
SPEOG ( I-J ) 

DH 


LENGTH+J2 
= SPYHFOG(M) 


PAGE 010 


14 


1 YMINf YMAXf 1 f 0) 


0407 

0408 

0409 104 

0410 

C 

0411 150 

0412 

C 

0414 

0415 

0416 

C 

0418 

0419 

0420 

C 

0422 

0423 


CAUL MPLOT ( ILFT +500 » ITOP+300f -1 ) 

TYPE 104 

FORMAT ('+ 15-PT FILTERED SLOW PHASE EQG DATA'" ) 
CALL COPY CO) 

IBLK * IBLK - 3 

IF ( IBLK .LT. NBLK ) GO TO 48 

TYPE *?' ANOTHER CHANNEL OF DATA? <YES=1>' 

ACCEPT * t JES 

IF (JES .EQ. 1) GO TO 42 

TYPE *»' TRY ANOTHER FILE? <1*YES> ? ' 

ACCEPT *, MORE 

IF (MORE .EQ. 1) GO TO 10 

STOP 

END 
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for F'rosram Unit FP ID 


Storaae msp 

Local Vsrublesp . PSECT $DA TA • Size 

Name T ype Offset Name Type 

ANRL R*4 000454 Env DENOM R*4 

DSKINC 1*2 000516 Eov DSUM R*4 

I'XSUH R*4 042470 DXYSUM R*4 

FPRMS R*4 042414 FORMS T R*4 

H R*4 042404 HHT R*4 

IRLK 1*2 042302 IBOT 1*2 

ILF 1*2 042532 ILFT 1*2 

IRIT 1*2 042344 IRMODE 1*2 

ISL 1*2 042510 ISLG 1*2 

ISTEP 1*2 042304 ISTP 1*2 

ITEMP 1*2 042310 ITLST 1*2 

IX 1*2 042540 IXS 1*2 

J 1*2 042330 JES 1*2 

K 1*2 042274 L 1*2 

LN 1*2 042460 LOREN 1*2 

M 1*2 042554 MORE 1*2 

NBLK 1*2 000530 £«v NCO 1*2 

NOUMMY 1*2 042276 NEH 1*2 

NURDS 1*2 042306 PERIOD R*4 

SLOPE R*4 042522 SSFD R*4 

SSHSD R*4 042440 SSSD R*4 

THT R*4 042266 XILF R*4 

XMAX R*4 042356 XMEAN R*4 

XMIN R*4 042352 XN R*4 

YINT R*4 042526 YMAX R*4 

Local and COMMON Arrays: 


= 042652 < £917. words; 


Off set 

Name 

T ype 

0 f f s e t 

042516 

DH 

ft* 4 

042550 

042544 

DXSSUft 

ft*4 

042474 

042504 

DYSUM 

ft*4 

042500 

042446 

FI2 

R*4 

042512 

0422 62 

I 

1*2 

042320 

042346 

I ERR 

1*2 

042250 

042342 

INT RAT 

ft*4 

0004 76 Ect v 

042312 

IS 

1*2 

042374 

042272 

ISN 

1*2 

042322 

042462 

1ST RT 

1*2 

042456 

042454 

I TOP 

1*2 

042350 

042542 

I YES 

1*2 

042362 

042444 

J2 

1*2 

042410 

042316 

LENGTH 

1*2 

042246 

042300 

LTERrt 

1*2 

042452 

042556 

N 

1*2 

0422c0 

042376 

NDIF 

Cl 

PH 

042412 

042256 

NPTS 

1*2 

042314 

000460 Eqv 

SDRHS 

R*4 

042430 

042420 

SSMFD 

R*4 

042424 

042434 

SUM 

R*4 

042400 

042534 

XLGT 

R*4 

042324 

042370 

XMFSD 

R*4 

042252 

042464 

XSUM 

R*4 

042364 

042336 

YMIN 

R*4 

042332 


Name Type 

BUFR2 1*2 
COEF R*4 
FDEOG R*4 
FILNAM L*1 
HE06 R*4 
ICHNUM 1*2 ’ 

IPOSN 1*2 
MAX 1*2 
SDEOG R*4 
SF'EOG R*4 
SPYHEO R*4 
XT R*4 

YHEOG R*4 
Z 1*2 


Subroutines » Functions r Statement and Processor-Defined Functions: 


Name 

Type 

Name 

Type 

Name 

Type 

Name 

Type 

Name 

T ype 

ARS 

R*4 

ANOTAT 

R*4 

CHRSIZ 

R*4 

COPY 

R*4 

DISKIO 

ft *4 

ERASE 

SORT 

R*4 

R*4 

FLOAT 

TSXCHK 

R*4 

R*4 

GRID 

XYF'LOT 

R*4 

R*4 

GRINIT 

ft*4 

MPLOT 

1*2 


Section 

Offset 


Si 

ze 

Dimensions 

♦DATA 

036076 

004000 

( 

1024 . ) 

(1024) 

♦ DATA 

012002 

000074 

( 

30 e ) 

(15) 

♦ DATA 

032076 

004000 

< 

1024 . ) 

(512) 

♦ DATA 

042076 

000014 

< 

6 « ) 

(12) 

♦ DATA 

002002 

004000 

( 

1 024 . ) 

(512) 

♦ DATA 

000542 

000040 

< 

16a ) 

(16) 

♦ DATA 

000702 

000040 

( 

16a ) 

(16) 

♦ DATA 

000742 

000040 

( 

16a ) 

(16) 

♦ DATA 

016076 

004000 

< 

1024a ) 

(512) 

♦ DATA 

022076 

004000 

( 

1024a ) 

(512) 

♦ DATA 

026076 

004000 

c 

1024a ) 

(512) 

♦ DATA 

006002 

004000 

( 

1024a ) 

(512) 

♦ DATA 

012076 

004000 

< 

1024a ) 

(512) 

♦ DATA 

000002 

002000 

( 

512. ) 

(512) 
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